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Abstract 

Neutrino telescopes will open, in the next years, new opportunities in observational 
high energy astrophysics. In these detectors, atmospheric muons from primary cos- 
mic ray interactions in the atmosphere play an important role, because they provide 
the most abundant source of events for calibration and test. On the other side, they 
represent the major background source. 

In this paper a fast Monte Carlo generator (called MUPAGE) of bundles of at- 
mospheric muons for underwater /ice neutrino telescopes is presented. MUPAGE is 
based on parametric formulas [1] obtained from a full Monte Carlo simulation of 
cosmic ray showers generating muons in bundle, which are propagated down to 5 
km w.e. It produces the event kinematics on the surface of a user-defined cylinder, 
surrounding the virtual detector. The multiplicity of the muons in the bundle, the 
muon lateral distribution and energy spectrum are simulated according to a specific 
model of primary cosmic ray flux, with constraints from measurements of the muon 
flux with underground experiments. As an example of application, the result of the 
generation of events on a cylindrical surface of ~ 1.4 km^ at a depth of 2450 m of 
water is presented. 

PACS: 95.55.Vj 

Key words: Simulation of atmospheric muons; neutrino telescopes; multiple 
muons; Monte Carlo generator. 



PROGRAM SUMMARY 



Title of program: MUPAGE 

Nature of the physics problem: Fast simulation of atmospheric muon bundles for 
^ Corresponding author 



Preprint submitted to Elsevier Science 



4 July 2008 



underwater /ice neutrino telescopes. 

Method of solution: Atmospheric muon events are generated according to parametric 
formulas [1] giving the flux, the multiplicity, the radial distribution and the energy 
spectrum. 

Licensing provisions: none. 
Programming language used: C++ 

Computers on which the program has been tested: Pentium M, 2.0 GHz; 2x Intel 
Xeon Quad Core, 2.33 GHz 

Number of processors used: one 

Operating systems under which the program has been tested: Scientific Linux 3.x; 
Scientific Linux 4.x; Slackware 12.0.0 
Number of bits in a word: 32 

Memory required to execute program with typical data: 50 MB 
Has the code been vectorized or parallelized? no 
Number of lines in distributed program: 1 933 
Number of bytes in distributed program: 59 166 
Distribution format: tar 
Library used: ROOT [26] 

Additional comments: The program requires the ROOT libraries for the pseudoran- 
dom number generator. 

Restrictions of the program: Water vertical depth range from 1.5 to 5 km w.e.; zenith 
angle range from to 85 degrees. 

Keywords: Simulation of atmospheric muons, neutrino telescopes, multiple muons, 
Monte Carlo generator 
PACS: 95.55.Vj 

Classification: 1.1 Cosmic Rays; 11.3 cascade and shower simulation. 
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1 Introduction 



Most astrophysical sources are expected to produce neutrinos in addition to 
photons and cosmic rays [2]. Theoretical predictions for neutrino fluxes indi- 
cate that ~ 1 km^ scale apparatus is probably needed. A neutrino telescope 
in the South Pole is currently under construction and is taking data with a 
reduced configuration [3]. In the Northern hemisphere, the European consor- 
tium KM3NeT [4], including the three collaborations, ANTARES [5], NEMO 
[6] and NESTOR [7] , is working on a design study for a large deep-sea neutrino 
telescope in the Mediterranean sea. 

Although neutrino telescopes are located under large depth of water or ice, 
a great number of high energy atmospheric muons can reach the active vol- 
ume of the detectors. Atmospheric muons are produced in the decay of charged 
mesons generated by the interactions of primary cosmic rays with atmospheric 
nuclei. They represent the most abundant signal in a neutrino telescope and 
they can be used to calibrate the detector and to check its response to the 
passage of charged particles. On the other side, they can constitute the major 
background source, mainly because downward-going muons can incorrectly be 
reconstructed as upward-going particles mimicking high energy neutrino inter- 
actions. Muons in bundles seem to be particularly dangerous [8]. A full Monte 
Carlo (MC) simulation, starting from the generation of atmospheric show- 
ers, can accurately reproduce the main features of muons reaching a neutrino 
telescope, but it requires a large amount of CPU time. 

In this paper an event generator (MUPAGE: MUon CEnerator from PAra- 
metric formulas) based on parametric formulas [1] is presented. The formulas 
allow to calculate the flux and angular distribution of underwater/ice muon 
bundles, taking into account the muon multiplicity and the multi-parameter 
dependent energy spectrum. The range of validity extends from 1.5 km to 
5.0 km of water vertical depth, and from 0° up to 85° for the zenith angle. 
MUPAGE output is an ASCII table containing the kinematics of events at 
the surface of a can. The can is an imaginary cylinder surrounding the active 
volume of a generic detector. The ASCII table can be used as input in the fol- 
lowing steps of a detector-dependent MC simulation, which include production 
of light in water/ice and simulation of the signal in the detection devices. 

This paper is organized as follows. In Section 2 the main features of the under- 
water muon flux, described by the parameterizations of [1], are reviewed. In 
Section 3 the MUPAGE structure is presented, in particular details are given 
on how to run the program, the input parameters and the output file format. 
Section 4 describes the generation of events on the surface of the can, with the 
use of a Hit-or-Miss method. Events can be single muons (Section 5) or the 
more complex multiple muons (Section 6). MUPAGE produces events with 
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the same relative weight. For each run the hvetime is computed as described 
in Section 7. Finally, as an example of application, the case of the ANTARES 
detector is presented in Section 8. 



2 Simulation of bundles of atmospheric muons 

Atmospheric muon flux in the deep water is at present experimentally poorly 
known [9]. Some parameterisations for the underwater flux and energy spec- 
trum are available in literature [10,11,12]. None of them gives the possibility 
to simulate two or more muons produced in the same cosmic ray interaction 
(muon bundles). To overcome this limitation, at a price of huge CPU time re- 
quirement, full simulations of atmospheric showers induced by primary cosmic 
ray interactions are performed by neutrino telescope collaborations [13,14]. 
These full MCs are generally based on the CORSIKA [15] or HEMAS [16] 
packages. They are then followed by the propagation of the surviving muons 
from sea level to the detector position. 

MUFAGE rehes on parametric formulas which describe the flux, the angular 
distribution and the energy spectrum of muon bundles for vertical depths h 
from 1.5 to 5 km w.c. of water or ice. The muon energy depends on h, on 
the zenith angle i) and, for muon in bundles, on bundle multiplicity m and on 
the distance R of each muon with respect to the shower axis. The parametric 
formulas were obtained starting from a full MC simulation of primary cosmic 
ray interactions and shower propagation in the atmosphere using the HEMAS 
code. HEMAS was preferred since it was largely used and cross-checked with 
MACRO data, a large underground experiment operating from 1994 to 2000 
and located at the INFN Gran Sasso laboratory, at a depth comparable to 
that of neutrino telescopes. In particular, MACRO measured the flux of muon 
bundles [17] and the distribution of distances between muon pairs in a muon 
bundle (the so-called decoherencc distribution [18]). The input primary cosmic 
ray spectrum used in MUFAGE was an unpublished MACRO model bounded 
by the measurements of underground muons. To optimize the full MC sim- 
ulation, the code was restricted to follow secondary particles with energies 
E > 500 GeV. Muons reaching the sea level with energies E > 500 GeV were 
propagated through water down to 5 km using the MUSIC (MUon Simulation 
Code) program [19]. The so called prompt muons (originated from the decay 
of charmed mesons and other short-lived particles) were not included. They 
were expected to give a not negligible contribution for muon energies from 
~ 10 TeV to ~ 10^ TeV, depending on the charm production model [20]. 

A MUFAGE event is a bundle of muons with multiplicity rric on the can. 
Muons in the bundle are assumed to be parallel to the shower axis, and to 
reach at the same time a plane perpendicular to the shower axis. The bundle 
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multiplicity, direction and impact point of the shower axis on the can surface 
are generated first. Then, for each muon in the bundle, the distance from the 
shower axis, the energy and the coordinates of the impact point on the can 
surface are calculated. 



2. 1 Parametric formulas for the flux, energy spectrum and radial distance of 
muons 



The flux of bundles with muon multiplicity m at a given vertical depth h and 
zenith angle 'd is parameterized using two free parameters {K and v) as: 



$(/i,^,m) = ^^%$ (1) 



As for the following formulas in this subsection, more details on the functional 
dependence of the parameters are reported in [1]. Fig. 1 shows the flux of 
muon bundles with different multiplicity at the depth h = 2.5 km w.e. for five 
different zenith angles, = 0°, 20°, 40°, 60° and 70° obtained using (1). 

Fig. 2 shows the flux of muons (with energy > 20 GeV) versus the zenith 
angle cosine at three different depths. In this case, a bundle of multiplicity 
m is accoimted for m muons. The full lines are obtained with (1), for depths 
h = 1.64, 2.0 and 3.0 km w.e. The points represent the AMANDA-II unfolded 
measurement [21], at the depth h — 1.64 km w.e. and with the same muon 
energy threshold. The (red) dotted hnes are from [12], at depths h — 1.61, 2.0 
and 3.0 km w.e. In [12], the bundle multiplicity was not taken into account 
and the depth-intensity relation of Bugaev et al. [11] was used as input. In our 
case, the depth-intensity relation for the vertical direction (called in literature 
l{h,0)) is computed with (1) and compared with other parameterizations in 
Fig. 1 of [1]. The (normalized) difference between what obtained with (1) and 
what reported in [11] is -16% at 1.5 km w.e. and -5% at 5.0 km w.e. The dif- 
ferences at the same depths of (1) with respect to the Okada parameterization 
[10] are +1% and -7%, respectively. 

The expected energy distribution of muons, assuming a power-law for the 
primary beam energy, at a slant depth X — h/cos'& is [22]: 

= G • E,e^''^'-^\E, + 6(1 - e-^^)]-^ (2) 



d{logwE^) 



The quantities e and 7 are considered as free parameters, while /3 as a constant. 
G — G(7, e) represents a normalization factor, in order to get the integrated 
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Fig. 1. Flux of muon bundles with different multiplicity m at the depth of h = 2.5 
km w.e. for five different zenith angles, 'd = 0°, 20°, 40°, 60° and 70° as computed 
with (1). 

muon energy spectrum up to 5 x 10^ GeV equal to 1. MUPAGE uses (2) to 
extract the muon energy, both for single muon events and for muon bundles. 

• Single muons (bundles with multiplicity m = 1). The parameter 7 
depends on the vertical depth h only: 7 = 7(/i). The parameter e depends 
both on h and on the zenith angle ■}}:€ = e{h,d). See Section 5. 

• Multiple muons (bundles with multiplicity m > 1). Muons produced 
in the decay of secondary mesons follow the energy distribution of the parent 
mesons. As a consequence, the most energetic muons are expected to be 
closer to the axis shower. To sample the energy of a muon in a bundle, its 
radial distance R from the shower axis is taken into account. The muon 
radial distance distribution is described as: 

^ = C ^ (3) 
dR [R + RoY 

where the parameter Rq depends on the vertical depth on the bundle 
multiphcity m and on the zenith angle Rq = Ro{h,m,{}); a depends on 
h and m: a = a{h,m). Finally, C = C{Ro,a) is a normalization factor. 
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Fig. 2. Muon flux (with energy > 20 GeV) as a function of the cosine of the zenith 
angle and for three different depths. The full lines are obtained from Eq. (1), for 
depths h = 1.64, 2.0 and 3.0 km w.e. from top to bottom. The points represent the 
AMANDA-II measurement [21], at the depths h = 1.64 km w.e. and with the same 
muon energy threshold. The (red) dotted lines are from [12], at depth h = 1.61,2.0 
and 3.0 km w.e. 

For each muon in the bundle, the energy is extracted with a probability 
given by the distribution (2). In this case, the parameter 7 depends on the 
vertical depth h, on the bundle multiplicity m and on the radial distance 
R: •y = 'y{h,m,R). The parameter e depends on h, on the zenith angle -d, 
and on R: e = e{h, {}, R). See Section 6. 

Eqs. (1), (2) and (3) are implemented in the MUFAGE Muons class as described 
in the next section. The values of the numerical constants are implemented in 
the Parameters class. 



3 Program structure 

MUFAGE code is written in C++ and it has been tested with gcc version 
3.2.x. , 3.4.x and 4.1.x. The program requires ROOT libraries [26] for the pseu- 
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C- Shell script file 


template script (tcsh) 


run-mupage.sh 


Shell script file 


template script (bash) 


rtJI/AUlVIJi/ 


AoUll me 


instructions on running MUPAGE 


dat 


directory 


input parameter files (see Section 3.1) 


evt 


directory 


event output files (see Section 3.2) 


inc 


directory 


C++ include files for MUPAGE 


livetime 


directory 


run livetime (see Section 7) 


src 


directory 


MUPAGE C++ source code 


Table 1 



Description of the files contained in the tar file main directory mupage/ 

dorandom number generator (see Section 4.1). Files from the tar archive are 
extracted in the main folder mupage/. This folder contains the Makefile, a 
README ASCII file, two template script files (for tcsh and for bash) to launch 
the executable and some subdirectories as described in Table 1. After running 
Makefile, a Linux/ directory is created, containing the object files; the exe- 
cutable is created in the main directory. The prototypes files ( . hh) are in the 
inc/ directory and the declaration files ( . cc) in the src/ directory. They do 
not need to be changed by the user: 

inc/ConvertingUnits .hh : units conversion 

inc/Geometry . hh : geometry definition 

inc/Parameters .hh & src/Parameters . cc: class with parameters from [1] to 

compute the parametric formulas 

inc/Decode .hh & src/Decode . cc : class to decode input parameters 

inc/Muons.hh & src/Muons . cc : class to generate an event 

src/mupage . cc : main program 



The simplest way to execute MUPAGE is the use of the C-shell script file 
run-mupage . csh for tcsh, or the Shell script file run-mupage . sh for bash. 
In both template scripts the user can modify random seed, run number and 
the number Ngen of events to be generated. It is also possible to define the 
directories for MUPAGE output files. The livetime of the run is evaluated 
from Ngen (see Section 7). The following control cards, defined in the template 
scripts, can be set as -character [name of data] : 
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-n 




help 


-i 


Jprun_ia 


run number 


-n 


$num_events 


number of events to be generated 


-s 


$seed 


random seed (default 0) 


-P 


SINFILE 


full name (including path) of the input file (see Section 


-o 


SOUTFILEl 


full name (including path) of the first output file 




$OUTFILE2 


full name (including path) of the second output file 
(for both output files see Section 3.2) 



Card values are provided guments to the executable: 

./mupage.exe -i $run_id -n $nuiii_events -s $seed -p $INFILE 
-o SOUTFILEl $0UTFILE2 



3. 1 Description of the input file 

The generator needs some parameters as input, which are contained in the 
file parajneters.dat in the directory dat/. Parameters refer to the detector 
configuration (see Fig. 3), to the medium and to the range of simulation pa- 
rameters (multiplicity, zenith angle, muon energy). Name, default value and 
unit of parameters are described in Table 

Note that, in the input file, the MULTmax parameter has a default value 
equal to 1000. This is the only parameter affecting CPU time significantly. To 
optimize CPU time, the event files can be produced with different ranges of 
[MULTmin, MULTmax]. In this case, when filling histograms, results must 
be normalized taking into account different values of livetime for each file (see 
Section 8) . 

3.2 Description of the output files 

The code provides two output data files: SOUTFILEl e $OUTFILE2. The 
former is written (by default) in the directory mupage/evt and contains all 
information about the generated events in a formatted ASCII table. An ex- 
ample of output data file with 1000 generated events (run_01 .evt) is present 

^ Default values refer to the ANTARES experiment in the Mediterranean sea, see 
Section 8. 
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Fig. 3. Sketch of the geometrical meaning of some input parameters. The cyhnder 
surrounding the instrumented volume is the can, with radius Rcan and height Hcan- 
The events are generated on an extended can, with R^xt = Rcan + Recr- The origin 
of coordinate system lies on the cylinder axis, but not necessarily at the center of 
the cylinder. The lower disk is at a depth Hmax with respect to the sea/ice surface. 
The origin of the detector coordinate system lies on the cylinder axis, but it does 
not necessarily coincide with the center of the cylinder. 



in the folder. Each line of the table contains the following information: 

eventJd mult track Jd Xi Ui Zi Vy Vz Ei ti fijid 

where: 

• eventJd = event number; 

• mult = multiplicity of the muon bundle at the depth where the shower axis 
hits the can; 

• for each muon in the bundle and intercepting the can {m^ < mult): 

- track Jd = i [i = l,mc), muon identifier in the event; 

- (xj, Hi, Zi), coordinates of the muon impact point on the can surface; 

- {vx, Vy, Vz), direction cosines of the muon, coincident with those of the 
bundle axis; 

- Ei (in GeV), energy of the muon; 

- ti, time delay of the i-th muon at the can surface with respect to the first 
muon in the list, {i = 1). ti can be either positive or negative; 

- fiJd, muon particle identification number, for the detector-dependent simu- 
lations following the MUPAGE step. By default it is inserted the GEANT3 
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MULTmax {nimax) 


1000 


mcLximum muon multiplicity 


GEANTid 


6 


muon ID (default: the GEANT3 id [23]) 


MFactor 


1.0 


multiplicative factor for special geometries 



Table 2 



Description and default values of parameters in the input file 
/dat/parameters . dat. 

code fi~=6. GEANT4 uses the Particle Data Group encoding and /x~=13 
[23]. By inserting the appropriate ID number for the muon, any other MC 
can be used. 

A second output file $OUTFILE2 is written (by default) in the subdirectory 
mupage/livetime and contains the livetime (units: seconds and days) for the 
simulated run. The livetime is given with its computed statistical error. The 
evahiation of livetime is described in Section 7. All events have the same weight 
(=1). The flowchart of the Monte Carlo program is shown in Fig. 4; details 
are described in Sections 4, 5 and 6. 
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=<i>(m,h,e) •S(e)-Ar2(e) 

(sez. 4 2, eq.8) 




Fig. 4. The flowchart of the MUPAGE event generator. The smooth-angle rectan- 
gles indicate the extraction of uniformly distributed random values. The decisional 
rhombuses in bold select values according to formulas reported in Section 2, with a 
Hit-or-Miss method. The procedure is iterated for Ngen events. 

4 Generation of muon bundles on the can surface 



4-1 Sampling the bundle direction and impact point on the can surface 



Muon bundles at ~ 2 km water equivalent depth can have multiplicity as large 
as 10^, and muons can be hundreds of meters far from the axis shower. On 
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the other hand, muons traveUing few absorption lengths far from the detector 
have a small probability to give a signal on photomultipliers. In order to accept 
peripherals muons in large bundles, the radius of the generation volume is 
increased by a quantity EnlargedCANr (Recr): specified in data cards. The 
radius of the generation cylinder is Rext — Rcan + Rear- If Rcan coincides 
with the radius of the cylinder surrounding the instrumented volume of the 
detector, it is recommended to define Recr — ^OXabs- 

As a first step, a generic bundle with muon multiplicity m* e ['mmin,'>TT'max], 
random zenith angle [(^miny^max] and azimuth angle 0*G [0°,360°] in the 
detector frame is generated. The values that are extracted from uniform distri- 
butions, as m, -i? and (f) at this step, are denoted with a *. The pseudorandom 
number generator used in the program is the Mersenne Twister algorithm [25] 
and it is included in the ROOT hbraries (TRandomS class) [26] . The axis of 
the bundle with intercepts the extended can in a random point of 

coordinates {x*, y*, z*), which are computed in the following way. The cylin- 
der projected area seen by the bundle is shown in Fig. 5. The impact point 
{Xr, Yji) is a random point on this plane. Only downward going particles are 
generated. It means that points with z* — Z^i^, on the lower disk of the 
can are not considered. and Ly in Fig. 5 are defined as — 2Rext and 
Ly = i^can sin -I- 2RextC0S'&*. The coordinates Xr and Yr can assume the 
values: 



— Rext < Xr < Rext (4) 
sin.* + «„.cos.*) < < ^sin.* + «„,cos.* (5) 



The point {Xr, Yr) on the plane perpendicular to the shower direction can be 
on the upper disk or on the lateral surface of the can. It lies on the upper disk 
(grey area of Fig. 5) if: 

YR>^sm^*-RextCOS^* (6) 



and 



Xl [YR-{H,,j2)sm^*f 

Rlxt {RextCOsd*? - ^' 



Returning in 3-D, the point coordinates are: {x*, y*, z*), with x* = —Xr, 

y* = - ^ tan^* and z* = 

cos -d* 2 

If Eqs. (6), (7) are not both true, the shower axis hits the lateral surface. 
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Fig. 5. Left: sketch of the plane 11 perpendicular to the shower axis. The muon 
bundle has zenith angle "& with respect to the detector z-axis. The interception 
point of the shower axis is uniformly distributed on the cylinder projection on the 
n plane. They will be generated outside the black region; the events in the grey 
area lie in the upper disk of the cylinder, while the remaining in the lateral surface. 

In this case, the intersection point has coordinates x* = Rext cos (p' , y* = 



Rext sm 0', z* = \ where 0' = 0* + -tt - 

^ sm-iy^ 2 / 

arccos ( — ). The points (x*, y*, z*) are distributed uniformly on the can 

V Rext / 

surface (with the exclusion of the lower disk). 



4-2 Hit-or-Miss method to sample the impact point 



The flux m), Eq. (1), decreases with increasing depth h, zenith angle "i? 

and muon multiplicity m, as shown in Figs. 1 and 2. The procedure described 
in Section 4.1 extracts uniformly h*, and m*. A Hit-or-Miss method [24] 
is used to reproduce the correct dependence of the number of events on these 
variables. For each set of parameters {h*, ■^*, m*), the number of events arriv- 
ing on the projected area S{'d*) in a small solid angle Afi('(9*) centred around 
■d* is computed: 

Nprojih*, m*) = <^ih*, m*) ■ S{^*) ■ AQ{^*) (8) 
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where 



(9) 



An{^*) = 27r[cos(i?* - 0.5°) - cos(i?* + 0.5°)] 



(10) 



A random number u is then generated, with: 



0<u< N, 



max — 



(11) 



N^ax corresponds to the set of values {h, i?, m) for which the function ^{h, -d, m) 
S{'d) ■ Afi('j?) is maximum. This function has a maximum for the minimum 

value of the detector depth {h = Hmin), corresponding to the can upper disk, 
and for the minimum value of the range of muon multiplicities (m = rrimin)- 
The maximization in terms of the ?9 variable is more complex, due to the not 
trivial dependence of m) • S{-d) ■ Afi(??) on 

In order to save CPU time, the maximum of ^{Hmin, nT-min) ■ S{'d) ■ Afl('d) 
is computed as the product of the maximum of the functions ^{Hmin, ^min) 
and S{d) ■ AVL{d). The former has a maximum in correspondence of 'dmin- The 



Using this approximation N^ax is evaluated as in (11). The parameter set (/i*, 
'd* ,m*) is accepted if: 



In the following, to simplify the notation, (/i, d ,m) will be used and the impact 
point coordinates become {x, y, z). 

If the can height is much larger than the disk diameter, Hcan ^ Rexti the 
approximation used for Nj^ax in (11) is not valid. It must be stressed that 
this is NOT the case for present neutrino telescope configurations. However, 
to correct the procedure, a multiplicative factor in the input parameters is 
introduced (default value =1). When the error message appears: 

ERROR! Nmax must be larger than Nproj 

and the program stops, the user must change MFactor in the input data file 
(parameters.dat) to a value larger than 1 (usually MFactor ~ 1.2 is enough). 
As an alternative, the value of R^cr can be increased. 



latter has a maximum, from (9) and (10), for zenith angle d' = arctan 




u < Nproj {h*,^*,m*) 



(12) 
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5 Single muons 



The underwater/ice flux of atmospheric muons is dominated (Fig. 1) by events 
reaching the detector with multiphcity m — 1, the so-caUed single muons. In 
this case the muon direction is assumed coincident with the shower axis and 
the impact point is {x, y, z). The arrival time of the muon on the can surface 
is t = 0. The muon energy E is extracted according to (2), whose parameters 
depend on the vertical depth h of the impact point on the can surface and on 
the muon zenith angle "d. 

A value of log;^Q E* is generated randomly between log^o Eyyim and log^^Q E„iax 
{Emin, Emax are given in the data cards). The value E* is accepted (or rejected) 
according to the Hit-or-Miss method: a random number u' is generated be- 
tween and the maximum of (2) at depth h and zenith 



The maximum of (2) occurs in correspondence of E^"'^ — ^^^^^^jy^- The value 
£■* is accepted if: 

u' <—^^(K^-E*) (14) 



6 Multiple muons 



6. 1 The radial distance of muons with respect to the bundle axis 



For events with muon multiplicity m > 1, the distance R of each muon from 
the bundle axis (in a plane perpendicular to the axis) is calculated, according 
to the radial distribution (3). R depends on depth /i, on bundle multiplicity 
m and on the zenith angle d. It is useful to define a new reference frame 
{Bundle Axis Frame, BAFrame) , where the zbaf-s^-^^s, coincides with the axis 
shower, see Fig. 6. Each muon is located in a point (X, Y) of the plane 11 
perpendicular to the axis shower. The distance of the point {X, Y) from the 
origin of the BAFrame is called Ri (the muon radial distance). R,* is sampled 
randomly between Rmin and Rmax (both values from data cards). The Hit-or- 
Miss method is used to accept (or reject) the value Ri*. A random number u" 
is generated between and the maximum of the lateral distribution function 
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(dN/dR) at the given h, '& and m. Ri* is accepted if: 

«"<^(/i,^,m,i?*) (15) 



The coordinates in the BAFmme are computed from the selected Ri as X ~ 
Ri ■ cos P and Y — R^- sin where /? is a random number between and 27r. 



6.2 Coordinates of the multiple muons on the can surface 



The shower axis intercepts the can in the impact point (computed in Section 
4.2) with coordinates {x, y, z)—{xsA,ysA,zsA)- Then, referring to Fig. 6, for 
each muon in the bundle: 

- (X, Y) = coordinates of the muon in the BAFrame; 

- {Xfj,, y^, z^) = coordinates of the muon in the laboratory frame; 

- {vx,Vy,Vz) — (sin 1? COS 0, sin i? sin 0, cos'j?) = direction of the muon in the 
laboratory frame; 

- {xi,yi,Zi) = projection of the point (a;^,y^,^^) along the shower direction 
on the can surface. 

When the point [X, Y) in the BAFrame is known, the point (x^, y^, z^) in the 
laboratory frame can be computed using a general matrix A resulting from 
the composition of three rotations [27]. In the so-called 'X-convention' the 
rotations are defined by the Euler angles ($,6,^), where the first rotation 
is by an angle $ around the z-axis, the second one is by an angle © e [0, tt] 
around the x-axis and the third one is by an angle \& around the z-axis (again). 
There is a univocal relationship between the three Euler angles and the zenith 
(t?) and azimuth (0) angles: 

$ = -7r/2, e = and * = + n/2. 

The transformation of the point {X, Y) in the BAFrame into the point {x^, y^, z^) 
in the laboratory frame is defined as: 











= A 


Y 









The coordinates of the impact point of each muon on the can are obtained 
using the projection of each point {xi^,y^,Zi^) along the direction {vx,Vy,Vz). 
This is done using the straight line defined by the three parametric equations 
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{xi, Hi, Zi) = (x^, y^, Zf,) + k{v^, Vy, v^), where k = {Z^ax - z^)/vz. The impact 
point of the i-th muon in the bundle is on the upper disk of the can if the 
straight hne intercepts the plane z = Z^^ax with a;^ + < i^ga-j. In this case, 
the coordinates are: 

(^i) Vi-i ^i) — i^n + kVx, Hn + kVy, Zjnax) (17) 

If x'f + yf > /?ext> impact point (xj, yi,Zi) lies on the can lateral surface or 
it does not intercept the can at all. 

The intersection of the straight hne with the lateral surface of the can (defined 
by equation x'^ + y'^ — B^^t) gives a second degree equation als? + 26A + c = 0, 
with a — v"^ + Vy, h — VxXfj, + Vyy^ and c = + — B!^xf The solutions are 

A± = & ± a/A ^ A = 6^ — ac. If A < the i-th muon does not intercept 
a 

the can. For A > the two possible impact points are: 

{xi, Zi) = (x^ + A+Vx, y,j, + A^+Vy, Zf, + A+vJ (18) 
{xi, yi, Zi) = (x^ + A_t;^, + A-t;^, + A^v^) (19) 

As atmospheric muons are downward going, the solution with the larger value 
Zi is chosen. If Zi < Zmin the i-th muon does not intercept the can. The number 
of muons intercepting the can surface determines the bundle multiplicity rric < 
m at the can. 



6.3 Arrival time of the muons in the bundle 



All muons in the bundle are assumed to arrive at the same time on the plane 
n perpendicular to the shower axis. In general, each muon reaches the can 
surface at a different time. The distance between the impact point of the i-th 
muon Pi{xi,yi, Zi) and the coordinates of that muon on the plane 11 in the 
laboratory frame P^,{x^,y^,, z^) is: 

d{Pi, P^) = sign ■ ^{xi - x^)^ + {yi - y^)^ + {zi - z^Y (20) 

The arrival time of the first muon in the list (i = 1) on the can surface is taken 
as = 0. All the remaining muons, labelled with i = 2, ...,mc, can intercept 
the can earlier (ti < 0) or later (tj > 0). The relative time is computed from the 
distance d{Pi,P^) defined by (20). A distance is an always positive quantity, 
but the evaluation of the relative delay between muons in the bundle requires 
the definition of the sign in (20). Referring to Fig. 6, if < Zi the distance is 
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Fig. 6. Lateral view of the extended can (the rectangle). The shower axis has direc- 
tion {vx,Vy,Vz) and intercepts the detector in the point (xsAiVSA, zsa)- rn muons 
are generated on the plane 11 perpendicular to the shower axis and nic < m inter- 
cept the can. A muon has coordinates {X,Y) in the BAFrame and {x^,y^, Zfj_) in 
the laboratory frame. The point {xi,yi, Zi) is the projection of the point (x^, y^, z^) 
along the /i direction on the can. The time delay of each muon is evaluated from 
the distance between the points {xi,yi,Zi) and {x^,y^, Zfj_). 

assumed positive {sign = 1), otherwise it is negative {sign = — 1). The delay 
ti of the i-th muon with respect to the first one is: 

^ d{Pi,P,)-d{Pi,P,i) ^21) 
* c 

where c = 2.99 x 10^ m/s is the light velocity. Since the distances can be either 
positive or negative, also (in ns) will assume either positive or negative 
values. 
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6-4 Muon energy for multimuon events 

The last step is the choice of the energy of each muon in a multimuon bundle 
according to the energy distribution (3). The muon energy extracted from 
this distribution depends on vertical depth, zenith angle, on the multiplicity 
of the shower and on the radial distance of the muon from the shower axis. 
The steps described in Section 5 for the energy of single muons are repeated 
for the evaluation of the energy of each muon in the bundle. 



7 Livetime of the simulation 

MUPAGE computes automatically the detector livetime corresponding to the 
number of generated events Ng^n on the can surface. The number of simu- 
lated events N{AVli) in a small solid angle Afij = 27r(cos't?ii — cos^92j), with 
multiplicity m = mmin and with shower axis intercepting the can upper disk, 
are evaluated in 33 bin© The expected rate of muon events with multiplicity 
nirain ou the can upper disk with area vri^g^.^ at the depth H^^n, and in the 
solid angle AQi is: 

iVMc(Aa) = ^Hmin, ^i, m^in) ■ S ■ AQ, {s-') (22) 

where = {'du + '(92i)/2, and S = nR'^^^cos'di (m^) is the projected area of 
the upper disk. The equivalent livetime for each bin is: 

T(Aa) = iV(Aa)/iVMc(Aa) (s) (23) 

The livetime with its statistical error is computed as the weighted average of 
the 33 values of T(Af2j) (which have the same value, within statistical errors), 
and written in the $OUTFILE2 file. 



8 An example of application: the case of ANTARES 

The ANTARES (Astronomy with a Neutrino Telescope and Abyss environ- 
mental RESearch) collaboration is operating the largest neutrino telescope in 
the Northern hemisphere in a site 2475 m deep, 40 km off La-Seyne-sur-Mer 

^ The first bin is 0° < < 10°, then 30 bins of 2 degrees are considered. The last 
two bins are 70° < "d < 76° and 76° < ?? < 85°, respectively. The bin size was chosen 
in order to have a constant or at least an adequate statistical sample in each bin. 
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(Prance). The detector (completed on May 2008) consists of an array of 12 
lines separated one from each other on the seabed by 60-80 m. The instru- 
mented area is ~ 0.06 km^ [5]. The simulation of atmospheric muons is one 
of the main task for ANTARES as for other neutrino telescopes. A full MC 
simulation is used, at the cost price of a large CPU consuming time. The 
ANTARES software tools are described in [13]. 

MUPAGE was used to produce a large sample of atmospheric muons. A data 
set with a livetime equivalent to one month of real data was generated with 
input parameters reported in Table 2. 358 files were created, each one smaller 
than 2 GB, with 10^ events/file and corresponding to 7260 ± 4 s. The CPU 
time required to produce a file (on a 2xlntel Xeon Quad Core, 2.33 GHz) 
was about 1 hour. It is a relatively large amount of CPU, but it is small when 
compared to other steps of simulation, namely the tracking and the Cherenkov 
light emission, which need ~10 times more CPU. Each file, after triggering 
and data conversion to a format equivalent to raw data, is equivalent to one 
run of 2.02 h of real data and it is used to test the detector response. The 
rate of generated events for a detector with parameters equal to the default 
value of Table 2 is 1240 Hz. This rate holds for muons above the threshold 
energy (20 GeV) at the surface of the generation cylinder of 1.4 km^ and it 
is much larger (a factor almost 1000) than the actually observed event rate 
after applying a realistic trigger algorithm. To be triggered, muons must pass 
sufficiently close to the instrumented detector volume to produce hits on a 
minimum number of optical sensors. Analysis of the comparisons between the 
MUPAGE, the full MC simulation and the data is in progress [29]. 

A much larger data set of atmospheric muons is needed for background study 
of the high energy neutrino {Ei, > 100 TeV) diffuse flux [14,28]. In this case, it 
is not necessary to simulate the bulk of lower energy muon bundles. A data set 
equivalent to one year was generated, using an optimized choice of multiplicity 
ranges, and a very conservative cut on the energy threshold Etht' = J2i=i E^^i > 
3 TeV (set by the parameter Ethreshold, see Table 2). With E^^r > 3 TeV, 
the rate of generated events on the default detector cylinder reduces to 27 Hz. 
The event multiplicity was divided in 6 sub ranges: m — 1 (4.7 Hz); m = 2 
(4.4 Hz); m = 3 (3.1 Hz); m = 4 ^ 10 (10.4 Hz); m = 11 ^ 100 (4.2 Hz) and 
m = 101 ^ 1000 (0.04 Hz). Each generated file with m = 1, 2, 3, [4 ^ 10] and 
[11 -j- 100] contains 10^ events (in order to have size < 2 GB). It respectively 
corresponding to 59, 63, 90, 29 and 62 hours of livetime. 5 x 10^ events/file were 
generated for m = [101 1000] corresponding to a livetime of 350 hours/file. 
The total number of generated files is 852. The CPU time required to generate 
each file ranges from 10 to 30 minutes using the same processor quoted above 
(it increases with the value of the minimum multiplicity). The total CPU time 
needed for the MUPAGE simulation of this data set was 232 hours. 
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9 Conclusions 



In this paper, a fast generator (MUPAGE) of tlie kinematics of atmospheric 
muon bundles is presented. As input, it uses parametric formulas for the flux 
of single and multiple muons valid in the range 1.5 < h < 5.0 km w.e. and 
< 85°. The energy spectrum of single and multiple muons is also simulated, 
taking into account the dependence of the muon energy on the shower multi- 
plicity and on the distance of the muon from the shower axis. The generator 
represents an useful tool for underwater/ice neutrino telescope to produce 
large amount of simulated data. As an example, the generation rate of atmo- 
spheric muon bundles on a cylinder with area of 1.4 km^ laying at a depth of 
2475 m and with total energy larger that 3 TeV is 27 Hz. The 8.5 x 10^ events 
corresponding to one year of data, were produced with 232 hours of CPU on 
2.33 GHz processor. 
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